CONSTRUCTION AND SHARP CONSISTENCY ESTIMATES FOR 
ATOMISTIC/CONTINUUM COUPLING METHODS WITH GENERAL 
INTERFACES: A 2D MODEL PROBLEM 
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Abstract. We present a new variant of the geometry reconstruction approach for the for- 
mulation of atomistic/continuum couphng methods (a/c methods). For multi-body nearest- 
neighbour interactions on the 2D triangular lattice, we show that patch test consistent a/c 
methods can be constructed for arbitrary interface geometries. Moreover, we prove that all 
methods within this class are first-order consistent at the atomistic/continuum interface and 
second-order consistent in the interior of the continuum region. 



1. Introduction 

Atomistic/continuum coupling methods (a/c methods) are a class of coarse-graining tech- 
niques for the efficient simulation of atomistic systems with localized regions of interest inter- 
acting with long-range elastic effects that can be adequately described by a continuum model. 
We refer to [6], and references therein, for an introduction and discussion of applications. 

In the present work we are concerned with the construction and rigorous analysis of energy- 
based a/c methods in a 2D model problem. Our starting point is the geometry reconstruction 
approach proposed by Shimokawa et al fTT] and by E, Lu and Yang [3] for the construction of 
"consistent" a/c methods in 2D and 3D. We propose a new variant of that approach to define 
a modified site potential at the a/c interface, which has several free parameters. We then 
"fit" these parameters so that the resulting a/c hybrid energy satisfies an energy consistency 



condition and a force consistency condition (see (2.6) and (2.7) for the precise definition of 
these terms; in the terminology of quasicontinuum methods our hybrid energy is free of ghost 
forces) . 

Explicit constructions along these lines can be found in [T7j for pair potentials and in [3] 
for coupling a finite-range multi-body potential to a nearest-neighbour potential, for high- 
symmetry interfaces. Our focus in the present work is the coupling to a continuum model and 
interfaces with corners; both of these cases are only briefiy touched upon in [3]. 

In recent years there has been considerable activity in the numerical analysis literature on the 
classification and rigorous analysis of a/c methods (see PQEl El [IDl [12] and references therein). 
Much of this work has been restricted to one- dimensional problems; only very recently some 
progress has been made on the analysis of a/c methods in 2D and 3D [5l [9l ITT] . 

The first rigorous error estimates for the method proposed in [3] (together with a wider 
class of related methods), in more than one dimension, are presented in [9] for 2D finite range 
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multi-body interactions. The work [9j assumes the existence of an interface potential so that 
the resulting a/c energy satisfies certain energy and force consistency conditions (a variant 
of the patch test) and then established first-order consistency of the resulting a/c method in 
negative Sobolev norms. 

Several important questions remain open: 1. It is yet unclear whether constructions of the 
type proposed in [3l[T7] can be carried out for interfaces with corners. 2. The error estimates in 
[9] contain certain non-local terms that enforce unnatural assumptions (e.g., connectedness of 
the atomistic region). 3. Moreover, this nonlocality causes suboptimal error estimates; namely, 
it destroys the second-order consistency of the Cauchy-Born model (see, e.g., [HlHIin]), and 
an unnatural dependence of the interface width enters the error estimates. (Moreover, we note 
that the error estimates in [TT] for a different a/c method are only first-order as well.) 

The purpose of the present work is to investigate for a model problem whether these re- 
strictions are genuine, or of a technical nature. To that end we formulate an atomistic model 
on the 2D triangular lattice with nearest-neighbour multi-body interactions (effectively these 
are third neighbour interactions), and construct new a/c methods in the spirit of [3l [17]. We 
then prove that the resulting methods are all first-order consistent in the interface region and 
second-order consistent in the interior of the continuum region, which is the first generalisation 
of the optimal one-dimenional result [ini Theorem 3.1] to two dimensions. 

Although it may seem restrictive at first glance to consider only nearest-neighbour poten- 
tials, we note that this is in fact an important case to consider. For example, bond-angle 
potentials (which are included in our analysis) usually consider only angles between nearest- 
neighbour bonds. More generally, multi-body effects are usually restricted to very small inter- 
action neighbourhoods, while long-range effects are often only displayed in pair potentials (in 
particular, Lennard- Jones and Coulomb), which can be treated, for example, using Shapeev's 
method [HI [151 E] . 

2. Atomistic Continuum Coupling 

2.1. Atomistic model. We consider a nominally infinite crystal, but restrict admissible dis- 
placements to those with compact support. Thus we avoid any discussion of boundary condi- 
tions, which are unimportant for the purpose of this work. 

Let Qg denote a rotation through arclength 7r/3. As a reference configuration we choose the 
triangular lattice (see also Figure [T]): 

£ := AZ^, where A := (01,02), 

Oi := (1,0)"^, and aj := Q;^~^ai,j G Z. 

We will frequently use the following relationships between the vectors a-,: 

Oj+e = Oj, OLj+z = ~0'ji and aj_i + Oj+i = aj for all j G Z. 

For future reference we also define a := (aj)^^]^, and Fa := (Faj)^^i, for F G M^^^. 

Our choice of reference configuration is largely motivated by the fact that C possesses a 
canonical triangulation (see Figure [T| and ^2.2[ ), which will be convenient in our analysis. 

The set of displacements and deformations with compact support are given, respectively, by 

"^0 := {u : C ^ : u{x) 7^ for at most finitely many a; G £}, and 
%■= {y:C^m?:y-ide%]. 
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Figure 1. The 2D triangular lattice and its canonical triangulation. 

We remark that deformations are usually required to be at least invertible, but that we avoid 
this requirement by making simplifying assumptions on the interaction potential. 

A homogeneous deformation is a map : C ^ M^, y^{x) := Fa;, where F G M^^^. We note 
that ^ unless F = I. 

For a map f : £ — )■ M'^, G N, we define the forward finite difference operator 

Djv{x) := v{x + ttj) — v{x), X & C, j & Z, 

and we define the family of all nearest-neighbour finite differences as Dy{x) := {Djy{x))^^i. 

We assume that the atomistic interaction is described by a nearest-neighbour multi-body 
site energy potential V G C^(]R2x6)^ ^ith V{sl) = 0, so that the energy of a deformation y G '3^q 
is given by 

^a(y) ■.= Y.V{Dy{x)). 
The assumption V{a.) = guarantees that <^a(2/) is finite for all y G 'S^q. 

2.2. The Cauchy— Born approximation. For deformation fields y G W^'°°(]R^; M^), such 
that y — id has compact support, we define the Cauchy-Born energy functional 

My)-= [ W{dy)dx, where 1^(F) := ^V(Fa), 

W G C3(]R2''2;]R), is the Cauchy-Born stored energy function. The factor := \/3/2 is the 
volume of one primitive cell of C, that is, Vr(F) is the energy per unit volume of the lattice FC. 

If 2/ G ^ is a discrete deformation, then we define its Cauchy-Born energy through piecewise 
affine interpolation: The triangular lattice C has a canonical triangulation JT' into closed 
triangles depicted in Figure [l] Henceforth, we shall always identify a function v : C ^ M.'' with 
its Pi-interpolant, which belongs to W^'°^(]R^; M'^). For a discrete deformation y G we can 
then write the Cauchy-Born energy as 

<^c{y) = I W{dy) dx=Y, \T\W{dTy). (2.1) 

where we define dxy '■= dy{x)\x^T and note that |T| = VLq/2 for all triangles T & ^ . 
Note that W{\) = and hence Sc{y) is finite for all y G 

Alternatively, can be written in terms of site energies, which will be helpful for the 
definition of a/c methods. Each vertex x E C has six adjacent triangles, which we denote by 
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• Continuum Node 
■w- Interface Node 
o Atomistic Node 



Figure 2. Atomistic- interface-continuum domain decomposition. 



Tjjj := conv{x, x + aj, x + Oj+i}, j = 1, . . . , 6 (cf. Figure [s]). With this notation, 

= Yl y^'iDyix)), where V^iDyix)) ■-= ^11 ^(^^.,,1/). (2.2) 
xg£ i=i 

Note that 6 C^(]R^^^) is well-defined since dx^.y is determined by the finite differences 
Djy{x) and Dj+iy{x). 

2.3. A/c coupling via geometry reconstruction. Let A G C denote the set of all lattice 
sites for which we require full atomistic accuracy. We denote the set of interface lattice sites 

by 

X := \^x E C \ A \ X + aj E A ioT some j G {1, . . . , 6}}, 

and we denote the remaining lattice sites by C := £ \ (^ U X); cf. Figure [2] 
A general form for the constuction of a/c coupling energies is 

^.M = E nDyi^)) + E + E ^^(^^(^))' (2-3) 

x£A xex zec 

where V^,x E X, are the interface site potentials that define the method (the atomistic site 
potential and the continuum site potential are determined by the atomistic model). 

For example, if we choose = V, then we obtain the original quasicontinuum method [8J 
(the QCE method). It is well understood that the QCE method suffers from the occurance of 
ghost forces, which result in large modelling errors [H El |71 [TH [T6] . 

In the following we present a new variant of the geometry reconstruction approach [31 [T7] 
for constructing V^. We define the interface potential as 

V^{Dy{x)) := y(7^.Dy(x)), (2.4) 

where TZx is a geometry reconstruction operator of the general form 

6 

TlxDy{x) := {TixDjy{x)).^^, and TZxDjy{x) := y^^Cxj^iDiy{x). (2.5) 

1=1 

Here (C'a;,i,j)^,j=i, a; G X, are free parameters of the method that can be determined to improve 
the accuracy of the coupling scheme. 

We use the acronym "GR-AC method" (geometry reconstruction-based atomistic-to-continuum 



coupling method) to describe methods of the type (2.3) where the interface site potential is of 
the form ( |2^ . 
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We aim to determine parameters Cx,j,i such that the couphng energy satisfies the follow- 
ing conditions, which we label, respectively, local energy consistency and local force consistency: 



V^{?Si) = Viysi) VF G M'^', Va; G X, and (2.6) 

2x2 



^2x2 

/ac(x;yF)=0 VFGM'^^ Vxg£, (2.7) 
where fs,c{,x] y) is the force acting on the atom at site x, initially defined by 

/ac(a:;i/):=-^^GM^ for y G ^^o; 
dy{x) 

however, we immediately see that /ac involves only a sum over a finite set of lattice sites, and 
hence the formula can be extended to all maps y : C ^ M?. In particular, (2.7) is a well-posed 



condition. Taken together, we call (2.6) and (2.7) the patch test. A hybrid energy t^ac of the 



form (2.3) is called patch test consistent if it satisfies both conditions. 

In the remainder of the paper, we will determine choices of the parameters C^j^i for general 
a/c interface geometries that give patch test consistent coupling methods. Moreover, we will 
prove that for all parameter choices we determine, the resulting a/c method is first-order 
consistent at the interface and second-order consistent in the interior of the continuum region. 
This extends the optimal ID result in [TH]. 



Remark 2.1. 1. To obtain a method with improved complexity one should use a coarser 
finite element discretisation in the continuum region. It was seen in [H] that the coarsening 
step can be understood using standard finite element methodology, and hence we focus only 
on the modification of the model, and the resulting modelling errors. 

2. Realistic interaction potentials have singularities for colliding nuclei, i.e., for deformations 
that are not injective. Clearly, our assumption that V G C^(M^^^) contradicts this. It is 
conceptually easy to admit more general site potentials in our work, however, this would 
introduce additional technical steps that are of little relevance to the problems we wish to 
study. □ 

2.4. Additional assumptions and notation. We use | ■ | to denote the £^-norm on M", 
and the Frobenius norm on M"^™. Generic constants that are independent of the potential 
(and the constants defined in the following paragraphs) and the underlying deformations are 
denoted by c. Although it is possible in principle to trace all constants in our proofs, it would 
require additional non-trivial computations to optimize them. 

2.4.1. Properties ofV. We define notation for partial derivatives of V, for g G M^^^, as follows: 

d,V{g):=^^^eR', and 9„y(g) := G for ^, j G {1, . . . , 6}, 

and similarly, the third derivative di^^kVis) ^ M^^^^^, which we will never use explicitly. We 
will frequently also use the short-hand notation 

K,, := djV{Dy{x)), Vt,, := 9,r((ar?/)a), and V^f,, := 5,V^(Fa), 

as well as analogous notation for second derivatives and for the site potentials V^, V\ and for 



V^'^^ which is defined in (3.4) 
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Interpreting the second and third partial derivatives as muhi-hnear forms we define the 
global bounds 



M2 := y2 ^^P dijV{g)[hi,h2], and 

'■' \hi\ = \h2\=l 

6 

M3 := y2 sup sup dij^kV{g)[hi,h2,h]. 

\hiHh2Hhs\=l 

With this notation it is straightforward to show that 

6 

y \diV{g) - diV{h)\ < M2 max \gj - hA, for g, h e M^xe^ ^2.8) 

' j = l,...,6 

We also assume that V satisfies the point symmetry 

V{i~g,^s)%i) = Vig) Vg G M^^^ (2.9) 

The following identities are immediate consequences of this condition: 

diV{fai)= -di+3V{fai), for z = 1, . . . , 6, F G M^^^ (2.10) 
di,V{fs.) = 9,+3,,+3^(Fa), for z, J = 1, . . . , 6, F G M2x2_ (2.II) 



We will prove results on the class Y, of all site potentials that satisfy (2.9) 

y ■= {V e C3(M2^^) I V satisfies ^ }. 



We will frequently use the following shorthand notation for partial derivatives of V, when 
there is no ambiguity in their meaning: 

v., := d.ViDyix)), V,, := 9,\/(Fa), Vt, := Va^y,, 

and analogous symbols for other potentials that we will introduce throughout the text. 

2.4.2. Linear functionals. For ?/ G and m G "^o we denote the directional derivative of <fa by 

{5<^^{y),u) := lim . 

We call dSa^iy) the first variation of Sg, and understand it as an element of . We use 
analogous notation for other functionals. This paper is largely concerned with establishing 
bounds on the modelling error 5<C(?/) — 5Sa_c{y)- 

To obtain sharp error estimates in W^'^-like norms, one needs to bound modelling errors in 
negative Sobolev norms, or, in our case, discrete verions thereof. Let £ : ^0 be a linear 
functional, and let ^ + ^ = 1, 1 < < 00, then we define 

:= sup 

\\du\\^p,=l 
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X + aA 7^ 

\ J- .1 




a; + 05 X + ae 



XT,biXT,& Xt,6,Xt,1 XtA,Xt,5 



X + ai 





Xt,1,Xt,2 Xt,3,Xt,4 Xt,2,Xt,3 



Figure 3. Convention for the symbols Tj. and xtj- 



2.4.3. Notation for the lattice and the triangulation. C is the set of vertices of and we 
denote the set of edges of ^ by with edge midpoints m/, / G ^ . 

For each vertex x G £ and direction Oj, let Tj. j := conv{x, x + a^, x + aj+i} G .^5^, j = 1, . . . , 6 
(see Figure [3]). The edge (x,x + aj) is the intersection of the two elements T^j- and T^j-i- 
Moreover, let xy^- G £ be the unique lattice point so that both xtj,xtj- + flj G T (again, see 
Figure [3]) . 

2.4.4. Discrete regularity. To measure regularity or "smoothness" of discrete deformations 
1/ G we first define the symbols 



\D'^y{x)\ := max \DiDjy{x)\, and \D^y{x) 

jj = l,...,6 



max \DiDjDky{x)\., for x G £. 

i,j,fc=l,...,6 



With mild abuse of notation, we then define the norms 

WD'^yWu^A) ■= \\\D'^y\\\p(A), and \\D^y\\i>p(^j^) := \\\D^y\\\iv(^A), 
for any Ad C and y dW^. If the label ^ is omitted, then it is assumed that A = C 

3. Construction of the GR-AC Method 

In this section we carry out an explicit construction of the GR-AC method. Our results are 
variants of results in [3], however, since our ansatz is different from the one used in [3|, and 
since we wish to be precise about the equivalence of certain conditions, we provide details for 
all our proofs. 

We assume throughout the remainder of the paper that the reconstructed difference TZxDjy{x) 
may depend only on the original differences Dj^iy{x), Djy{x), and Dj+i?/(x), that is, 

Cx,j,i = for I (i — j) mod 6| > 1, i, j G {1, . . . , 6}, xEZ. (3.1) 



For future reference, we call (3.1) the one-sidedness condition. 



In ^3.1 and ^3.2 we derive general conditions on the parameters that are independent of the 



choice of the atomistic region. In ^3.3 and ^3.4 we then compute explicit sets of parameters. 
3.1. Conditions for local energy consistency. We first derive conditions for the local 



energy consistency condition (2.6). 



Proposition 3.1. Suppose that the parameters C^j^i satisfy the one-sidedness condition (3.1), 
then the interface potential satisfies the local energy consistency condition (2.6) for all 
potentials V if and only if 



C 



1-C 



for j 



.,6. 



(3.2) 
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Proof. We require that V^{Fsl) = l^(Fa), for arbitrary V E which is equivalent to 

6 

= X] for J = 1, • • • , 6. 



i=l 



b2x2 



and in view of (3.1), we obtain the condition 



Since this has to hold for arbitrary F G 

Since aj = aj_i + Oj+i, this is equivalent to 

(^a^jj-i + Cxjj — + {Cxjj+i + Cxjj — l)aj-i = 0, 

and since aj_i,aj+i are linearly independent, we obtain the condition that 



a 



Subtracting these two conditions gives Cxjj+i = C^jj-i, and hence we obtain (3.2). 



□ 



As a consequence of Assumption (3.1), and Proposition 3.1 we have reduced the number of 



free parameters to six for each site x G X. To simplify the subsequent notation, whenever the 

(3.3) 



parameters Cx,j,i are chosen to satisfy (3.2), we will write 



a 



Cx 



and note that C 



a 



i-c. 



x,3- 



Since it is equivalent to (2.6) we call (3.3) the local energy consistency condition as well. 



3.2. Conditions for local force consistency. We rewrite <fac in terms of a hybrid site 
potential 

r v\^), xgc, 

S,M = Y.'^:'^Dy{x)), where V^^^) := I Vi(g), x G X, (3.4) 

x&c y ^(g)) X E A. 

Lemma 3.2. Suppose that the parameters {Cx,j,i)^j=i,x G X, satisfy the one-sidedness 
condition (3.1) and local energy consistency (3.3). Moreover, let 

(3.5) 



1 for X E A and Cxj := 2/3 for x E C, 



1,...,6, 



and let {Cxj,i)fj=i, x E AUC, be defined to be compatible with (3.1) and (3.3); th 



len 



6 6 



f'^ix- Fid) = Y.^Cx-a,,,, - Cx,,)V^,j Vx G C. 
j=l i=l 



(3.6) 



Proof. Using the notation (3.4), we have 



(5^ac(Fid),w) = ^^9,V;^^(Fa) ■ DMx), 



xGC i=l 



which immediately gives 



- r(x; Fid) = J2 [^^Vx^-aS^^) - ^^V:'i^^ 



i=l 



(3.7) 
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• Continuum Node 
<r Interface Node 
o Atomistic Node 

^ o o o o o 
o o o o o 

Figure 4. The flat interface case. 



Witli tlie notation introduced in (3.5), we obtain 

6 6 



(Fa) ■ DMx) = J2 ^F,, Yl C.,,,D,u{x), 
j=i i=i 



1=1 



wliicli implies 



Combining (3.8) with (3.7) yields (3.6). 

Testing (3.6) for all V ^ Y and F G M^^^, we obtain the next result. 



□ 



Lemma 3.3. Suppose that the parameters {Cx,j,i)ij=i,x E I, satisfy one-sidedness (3.1) and 
local energy consistency (3.2). Then <fac satisfies local force consistency (2.7) for all V E Y if 
and only if 



^ ] {Cx-~ai,j,i ~ Cx-ai,j+3,i ~ CxJ,i + Ca;j+3,i) — V J — 1, 2, 3, Vx G C 



(3.9) 



i=l 



Proof. Using (3.6) and point symmetry (2.10) one readily checks that (3.9) is sufficient for 



force consistency (2.7). To show that (3.9) is also necessary we test (3.6) with 



which clearly belongs to the class Y, to obtain 

6 

-r (x; Fid) = Y Y.^Cx-a^,^, - a,.)(F - l)a, 

j=l,A i=l 
6 

= ^ ] iCx-ai,l,i ~ Cx-ai,i,i ~ Cx^l/i + Cx,4,i) (F — 1)^1. 
i=l 



For this expression to vanish for all F G M?^"^ we obtain precisely (3.9) for j = 1. For j 
the same argument applies. 



2,3 
□ 
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3.3. Explicit parameters for flat interfaces. We now give a characterisation, for a flat a/c 



interface, of all parameters satisfying the one-sidedness assumption (3.1), which give a patch 
test consistent a/c method. 



Proposition 3.4. Suppose that A = {x E C\x2 < 0}, X = {x G £ | X2 = 0} and C = {x E 

£ I a;2 > 0} (see Figure^. Then the parameters (C*z,j,i)fj=i; x El, satisfy the one-sidedness 



condition (|3.l|), energy consistency (3.3), and force consistency (3.9), if and only if 

CxA 



c. 



Cx 



+ai,4 



+ai ,j 



Vx eX, 

Vx GX, 



and 

J G {2,3,5,6}, 



(3.10) 
(3.11) 



where we have used the reduced parameters defined in (3.3). 



Proof. One-sidedness (3.1 ) and energy consistency (3.3 ) yields the reduced parameters {Cxj)j=i, 
X G X, satisfying (3.3). Recall also the extension (3.5) of these parameters for x G C. 

Let X+ := {x + a2 I X G X} and X_ := {x — 02 | x G X}. Clearly, we need to test (3.9) only 
for X G X U X_ U X+. Exploiting the symmetries of the problem it is also clear that we only 
need to consider j = 1,2. 

It is straightforward to verify through direct calculations that any set of coefficients satisfying 



Let j 



G X then we obtain that (3.10) is necessary from the force consistency 
condition (3.9), applied at x + 02 or x + a^. Let j = 2, then we obtain Cx,2 



(3.10), (3.11) satisfies the equivalent force consistency condition (3.9). 
and X 

2 



Cx 



x+ai ,2 



from 

the force consistency condition (3.9) applied at x + 02. Therefore, (3.10) and (3.11) are also 
necessary. □ 



Remark 3.1. We observe that the coefficients (Cx,j,j)fj=i, x G X, are not unique, but that 
we have considerable freedom in the construction of the GR-AC method: For each direction 
that is not aligned with the interface, there is a free parameter, while for each edge (x, x + ai) 
lying on the interface, there is one additional free parameter. This freedom will be reduced in 
the case of corners. □ 



3.4. Explicit parameters for general interfaces. For general interface geometries we make 
the following separation assumption. This assumption requires that, if the atomistic region 
can be decomposed into several connected components, then they must be separated by at 
least four "lattice hops". 



Assumption 3.5. Each vertex x G X has exactly two neighbours in X, and at least one 
neighbour in C. 

As in the fiat interface case, we can completely characterise all parameters within the one- 
sidedness assumption, which satisfy the patch test. 



Proposition 3.6. Let A C C be defined in such a way that the interface set X satisfies 
Assumption 3.5. and is not planar. Then the parameters {Cx,j,i)^j=i, x G X, satisfy the one- 



sidedness condition (3.1), energy consistency (3.3), and force consistency (3.9), if and only 
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if 



Cxj — C. 

Cx, 



x+aj j+3 



Va; G X, X + aj G X, 



C. 



x,3 



1 

2/3 



Va; gX, 
Vx gX, 



X + ttj G A, and 
X + G C, 



(3.12) 
(3.13) 
(3.14) 



where {Cx,j)j=i, x G X, are the reduced parameters defined in (3.3). 



Proof. As in the flat interface case, one-sidedness (3.1) and energy consistency (3.3) are equiv- 
alent to having the reduced parameters {Cxj)^j=i,x G X, satisfying (3.3). Recall also the 



extension (3.5) of these parameters for x G ^ U C. 

Let X+ := {x G C I 3aj, x + aj G X} and X_ := {x G ^ | 3aj, x + aj G X}. We need to test 
(3.9) only for x G XUX_ UX+. The necessity of (3.12) follows as in the fiat interface case. The 



necessity of (3.13) and (3.14) can be obtained by testing the corner sites in X± in the interface 



geometry depicted in Figure |2j 

To see that (3.12)-(3.14) are also sufficient one notes, first, that the corresponding coefficients 



always provide zero contribution on each edge for the sum in (3.9). Computing the force at 
X G X+ we see that the contribution from is the same as from V^, 
since the pure Cauchy-Born model passes (3.9). For x G X_ 



and must therefore cancel, 
the same argument applies. 
It remains to test (3.9) for x G X, at corners. Since (3.9) is a local condition, and due to 



Assumption |3.5[ one may assume that the interface has only one corner. Since all other sites 
are in equilibrium, and since the forces are conservative, it follows that the corner must also 
be in equilibrium. 



(Alternatively, one may check (3.9) through explicit computations for the corner geometry 
shown in Figure [21 All other geometries can be reduced to this one by symmetry.) □ 



Remark 3.2. We observe that, for a general interface, we only have freedom to choose the 
geometric reconstruction parameters along the interface, namely, for each interface edge there 
is one free parameter. □ 

4. Consistency of the Cauchy-Born Approximation 



Before we embark on the analysis of the GR-AC method (2.3), we establish a sharp consis- 
tency estimate for Cauchy-Born approximation. Related results were established in ^ , which 
require more stringent conditions on the smoothness of the deformation field. For the analysis 



of a/c methods a sharp consistency estimate, such as Theorem 4.2, is useful. In the remain- 
der of the section we establish technical results that are useful for the subsequent consistency 
analysis of the GR-AC method. 

4.1. Second-order consistency. A natural way to represent the first variation of S.^ is 

6 6 

(5^a(l/),w> = 5]5^9,V(Z}i/(x)) ■ D^u{x) = Y.Y.^x, ■ D,u{x), (4.1) 

x^C j=l xGC j=l 

where we use the notation Vxj := djV{Dy{x)). This representation can be interpreted as a 
sum over mesh edges. By contrast, the most natural representation of SS'c is 

(54(2/), = Yl \T\dW{dTy) : dru. (4.2) 
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To estimate S<^a_ — bS^ we will rewrite (4.2) in a form mimicking (4.1). The opposite approach 



is also possible, but does not lead as easily to second-order consistency estimates. 



Lemma 4.1. For y e %,T e ST , let Vtj := djVidry ■ a); then 

3 

(54(y),M> = 5^5^(\/t.,^,, + \/t.,^_„,) ■D,n(a;), W e %, and (4.3) 

xec j=i 

3 

M> =J2J2 (^-.J- - ^-+-.-^+3) • D^x), yu e %. (4.4) 

xeC j=l 



Proof. It is easy to see that 



1 ^ 

dW{F) = -J2d,V{F8i)0a„ 



and hence, using Qq = 2\T\ and dxu ■ aj = Dju{xT,j), 



^6 1 ^ 



Te5^ j=i 



(4.5) 



Every edge appears twice in this sum since it is shared between two elements; hence we obtain 
the edge representation 

6 

{5S^{y),u) = EE K^T.,. + ^T.,_,,) • D,u{x) Vn G (4.6) 



Since Dj+su{x + aj) = —Dju{x), and using Vtj+s = —Vtj (see ( 2.10[ )) we can reduce this 
sum as follows: 



x£C j=l 
3 

xi^C j=l 



This concludes the proof of (4.3) 



For the proof of (4.4) one only needs to use the identity Djj^^u{x + aj) = —Dju{x). □ 



Theorem 4.2. Let y G then 
where M2,M^ are defined in ^2.4-1 



(4.7) 
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X + ai 



Figure 5. Visualisation of the proof of Theorem 4.2 



Proof. It is useful to visualize this proof using Figure [5| and Figure |3] for additional detail. 
From Lemma 14.11 we obtain 

3 

{S^M-6^,{y),u) = ^^5,(x)-D,n(a;), (4.8) 

X&C j=l 

where 6j{x) := V^j - - Vt^^^j - VV,,,_ij. (4.9) 

In the following we estimate 6i{x) only; the remaining estimates follow by symmetry. 

Let F+ := dx^^y and F_ := dr^ ^y, then Vr^ -^^i = Vf+,i and Vt^ q^i = Vf_,i. Moreover we can 
Taylor expand 

6 

K,i = Vf+,1 + ^ VF+,ii(A2/(a;) - f+ai) + 0{\D^y{x)\'^), and similarly 

1=1 

6 

-K+ai,4 = - Vf_,4 - ^F_,4i(Al/(x + a,) - F_a,) + 0(1 Ay(x)|2) 

6 

= Vf_,i - J2 ^F_,i(.+3)(Al/(a; + ai) - F_a,) + 0(| A^(x)p) 

6 

= Vf_,i + 5^ rF_,H(-A+3y(x + ai) - F_a,) + 0(| A|/(x)|2). 



A careful analysis of the remainder shows that 0{\D^y{x)\'^) < ^J2ij=i l^'^vi^)? 
for some G M^^^. In the remainder of the proof we will suppress the argument 6. 

Clearly, Vp.^ii — Vf_^_^ii = 0{\D'^y{x)\) < X]j=i l^iy^l l-D^l/(a;)|, and hence we can deduce 
that 

6 

^li^) = YVF+,ii{Diy{x) - F+Oi - Di+?,y{x + oi) - F_ai) + 0[\D'^y{x)\^) 
1=1 

6 

= ^VF+,ii(Az/(x) - Al/(a;+) + Al/(x + ai - a,) - Al/(a;-)) +0(|A?/(x)|2) 

i=l 
6 

=: 5^VF„H5, + 0(|Ay(x)p), 
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where := xr^ -^^j and x~ := xr^ g^j. (These are simply the vertices in the two adjacent 
elements such that the identities f±ai = Diy{xf) hold.) 

Tracing the previous Taylor expansions, we see that, in the last estimate, 0{\D^y{x)\'^) < 
^Eh=i\du,V\ 

We compute £3 in detail but only give the results for the remaining coefficients: 

^3 = D-iy{x) - D^y{x + ai) + D:iy{x + ai - 03) - D3y(x - 03) 
= - DiD^yix) + DiD^yix - 03) = D^D^D^y^x). 
By performing similar calculations for i = 1, 2, 4, 5, 6, one finds 

ei = 82 = e& = 0, £4 = DiDiD^y^x), and £5 = DiD2D^y{x)] 

hence we obtain that 5j{x) = 0{\D'^y{x)\'^ + \D^y{x)\) (recall that we assumed, without loss 
of generality, that j = 1), where 0{\D^y{x)\) < J2i=3,4,5 \9i,iV\\D^y{x)\. 
Combining these estimates, we obtain 



3 i/p ^ ^ 



x£C j=l x£C j=l 

Elementary estimates yield 

3 i/p 



' < M2\\D^yh. + M,\\D^y\\l,, and 



xec j=i 
3 



(EEi^.-wr')* s (2v^)'"'(E ml*"!'')'"'. 

x&C j=l T£.T 

from which the result follows immediately. □ 



In the following subsections, we derive technical results related to Theorem 4.2, in prepara- 
tion for the proof of consistency of the GR-AC method. 



4.2. Stress tensors. We can re-interpret Theorem 4^ in terms of a second-order error esti- 
mate for certain stress tensors. If, for some y G there exist tensor fields Sa(?/; •), ^dv, •) £ 
Po(^)^^^, which satisfy the identities 



(5^a(2/),M) = J2 l^|Sa(2/;T) : dru, and (4.10) 
iy),u) = Y,\T\^ciy;T):dTU (4.11) 



re 5^ 

then we call Sa an atomistic stress tensor and a continuum stress tensors. 



It follows from (4.1) and (4.2) that 



1 ^ 

Ea(y; T) := — V V^^^j ® a„ and (4.12) 

1 ^ 

Eliy; T) := dW{dTy) = 7^ E ^^.^^ ® «i (^.13) 
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Figure 6. Notation for neighbouring triangles of T G S/' . 



satisfy (4.10) and (4.11), respectively. As we will see immediately, they are not the unique 
choices. 

In the following calculation (and later on as well) we denote by the unique neighbouring 
element of T e which shares an edge with direction with T; see Figure |6] 

With this notation, and using the fact that Dju{xT,j) = Dju{xTj,j), we observe that 



1 ^ 

(54(2/), u)= E 1^1^ E ■ DM^T,j) 

° j=l 

T&sr ° j=i 

Y.\T\{^il\{'^T, + VT^,)®aX:drU 
° 1=1 ^ 



(4.14) 



0; 



which yields the alternative continuum stress tensor 

6 



(4.15) 



Furthermore, if we write the Cauchy-Born energy in terms of the site energy (2.2), and 



apply the procedure used to derive Sa, then we obtain a third variant of the continuum stress 
tensor: 

6 



(4.16) 



We see that stress tensors are not uniquely defined by (4.11) and (4.10). This causes analjd}' 



ical difficulties when deriving consistency error estimates, which strongly depend on the choice 
of the stress tensors. For example we will show in the following result that is second-order 



consistent. By contrast, and are only first-order consistent (cf. Remark 4.1) 



(4.17) 



Lemma 4.3. Let y G then 

\j:,iy;T) -El{y;T)\ < c{Ms\D'y{x)\' + M2\DS{x)\) 
for allT e £^,x e T. 

Proof. This estimate is obtained by reversing the construction of E^ in (4.14), and applying 



the estimates obtained in the proof of Theorem 4.2 



□ 
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Remark 4-1- 



Taylor expansions show that S^, k = 1,3, are only first-order consistent, 
\E,{y-T) -J:':{y;T)\ < cM2\DMx)\ for x G T, 



but that a second-order estimate such as (4.17) would be false. The first-order estimate can 
also be obtained from the fact that Sa(?/F; •) = ^civf, •) = dW{F) for all F G M^''^. □ 

4.3. Divergence-free stress tensors. In the previous subsection, we have seen that the 



stress functions defined in (4.10) and (4.11 ) are not unique. It is therefore crucial to characterize 



all divergence-free tensors, which is the purpose of the present section. We call a piecewise 
constant tensor a G Po(^)^^^ divergence free, if it satisfies 



/ a : Swdx = ^ |r|a(r) : 9tM = Vm G '^c- 



(4.18) 



Ni(^) : = 



Divergence-free tensors can be characterised as 2D-curls of non-conforming Crouzeix-Raviart 
finite elements. Let Nl(=^^) be defined by 

'^lint(T) is linear for each T & ,^ 
V is continuous at all edge midpoints 

The degrees of freedom for functions w G Ni(^) are the nodal values at edge midpoints, w{qf), 
f G and the associated nodal basis functions are denoted by (f. 

We have the following characterization lemma [13] for divergence free tensor fields. Although 
we will never use the equivalence of the characterisation explicitly, it motivates much of our 
subsequent analysis. ***!*** 



Lemma 4.4. A tensor field a G Po(^)^^^ is divergence- free (i.e., satisfies {A.18)) if and 
only if there exists ip G Ni(^)^, such that a = difj], where J is the rotation by 7i/2. 

Proof. It is easy to show that every tensor of the form a = dwJ, w G Ni(^)^ satisfies (4.18), 
by checking the result for a single nodal basis function ip = Q. 

To show the reverse, let be a simply connected domain, which is a union of triangles 
T G Suppose that the number of vertices in VL is H^V , the number of interior vertices is 
#V/, the number of edges in VL is and the number of triangles in VL is t^T. 

We test ( |4.18 ) for all m G that are non-zero only in the interior of Vt. The dimension of 
all cr G Po(^^)^^^ satisfying (4.18) for those u can be at most A^T — 2^Vi. On the other hand. 



the dimension of Ni(f2)^ is 2^E and the dimension of rotated gradients of Crouzeix-Raviart 
functions, denoted by 9Ni(i7)^J, is 2^E — 2. We will show below that the following formula 
holds: 

4#T - 2#y, < 2# - 2, (4.19) 

which immediately implies that the subspace of divergence-free tensor coincides with (9Ni(f2)^J. 
Moreover, the representation is of course unique (up to a shift) and therefore independent of 
the choice of the domain. 



To prove (4.19), we use Euler's formula, 

#V - #E + #T 



(4.20) 



and the identify 

3#F = 2#E -W + Wi, (4.21) 

which is obtained by a simple counting argument. (Note that j^V — j^Vj is the number of 

boundary edges.) Subtracting (4.20) from (4.21) yields (4.19). □ 
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Continuum stress tensor correctors. We have different forms of continuum stress 
and S^, wfiicli all can be used to represent 5Sc in the form (4.11), and hence their 



differences must be divergence free. Lemma 4.4 characterises the form of these differences and 
motivates the following result. 

Lemma 4.5. Let y G then there exists a corrector ip'^^ly; •) G Ni(<^7)^ satisfying the 
following two properties: 

Corrector property: T.l{y; T) - ^Ky; T) = dip^^{y; T)J VT G ^; (4.22) 

Lipschitz property: m/) | < |M2||-D^|/||£oo(jn£) V/ G ^. (4.23) 



Proof. Property (4.22) follows of course from Lemma 4.4, however, to establish (4.23) we 
require an explicit expression of ip"^^. We give the details of the proof for the case of an upward 
pointing triangle T G ^ (cf. the left configuration in Figure [6|. An elementary computation. 



1 



starting from (4.15) and (|4.16|) and using the symmetry property (|2.10|), yields 

i:liy;T)-J:liy;T) 



[{Vt,i - Vt„i) + {Vt,3 - Vt„3) + {Vt,5 - Vt„5)] ®ai + 



where ". . . " stands for terms that are symmetric to the ones in the first line. The directions 
ai, 03, as are chosen anti-clockwise with respect to the element T. 

We now observe that, if / is an edge of T with direction aj, j G {1, 3, 5}, then 



no 1 ' 



in T, 



in Tj. 



Let / be the edge of T with direction oi, then choosing 

^p'^^{y■,mf) 



o o 



(4.24) 



(4.25) 



and making analogous choices for the remaining edges, we obtain (4.22). 

With this explicit representation we can now prove the Lipschitz property (4.23). Let / 
denote the edge of T with direction ai, F := dry and Fi := dr^y, then 



|^''(2/; mf)\ < l\Vf,,i - Vf,i\ + l\Vf,,2 - V,,2\ + l\V,„3 - Vf,3\ 

6 

< IE (l^yl + 1^2,1 + lK3,l)|(Fi - F)a,| < .maxJ(Fi - F)a 



(4.26) 



where Vij = dijV{Gj ■ a) for some Gj G M^^^, and M2 is defined in ^2.4.1 One now verifies 
that 

(Fi - F)ai = 0, (Fi - F)a2 = DQD2y{xT,i), and (Fi - F)a3 = D^D3y{xT,b), 
which implies 

max |(Fi — F)aj| < max i\D'^y{xT \D'^y{xT a)\) ■ 

j=l,...,6 



Combining this estimate with (4.26) we obtain (4.23) for edges aligned with Oi. The remaining 
cases follow from symmetry considerations. □ 
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5. Consistency of the GR-AC Method 



We are now ready to state the second main result of this paper. The proof is estabhshed in 
^ through §5.3[ For the remainder of this section we assume that the hypotheses stated in 
Theorem 15.11 hold. 



Theorem 5.1. Let S'^c be defined by (3.4), with parameters (C^,ij )f x G X, satisfying the 
one-sidedness condition (3.1), as well as the patch test conditions (2.6) and (2.7). Suppose in 



addition that the parameters are bounded, that is, 

sup max =. C < +oo. 

xex i.«e{iv,6} 

Then there exists a constant Cx = Cx{C), such that 

||5^ac(z/) - 5S^{y)\l^,-,., < c(CxM2||L''|/||,P(xc..) + M2\\D^yhpic) + M4D''y\\%,^c)), 
where I^^^ := {x E C \ dist(x,X) < 1} is an extended interface region. 



(5.1) 



5.1. An a/c stress tensor. Following the construction of Sa in (4.12) (with replaced by 
<^ac), we obtain a representation of 5<fac in terms of an a/c stress Sac: let y E and u E '^o, 
then 



[S^i,ciy),u) = ^ |T|Sac(l/;r) : dru, where 



1 6 



(5.2) 
(5.3) 



and we recall that V^^^ = djV^'^{x; Dy{x)). We now require the following additional notation: 



=^1 



{T G ^|Tn (XUC) = 0}, 
{Tg ^ITn(xu^) = 0}, 

Sr\{Src\^ STj^), and 



^\(^cU^^). 



(5.4) 



Lemma 5.2. (i) Let Sac be defined by ( |5.3 ), then, for all y E 

Sac(i/;T) = Sa(z/;r) VT G 5^1, and (5.5) 

Sac(i/;T) = S3(t/;T) VTGe^c- (5.6) 

C^j Lei F G M^^^,- i/ien t/iere exzsfo a unique ilj'^iV; •) G Ni(e^)^ such that 

Sac(z/F;T)-Sa(i/F;X) = 9V^^^(F;T)J VTg.^, and (5.7) 

^'^^(F; m/) = V/ G U ^c- (5.8) 

Moreover, there exists Lac depending only on C such that the following Lipschitz property holds: 

\r'{f; nif)- r'{G; mf)\< L^^M^IF - G| VF, G G M2^^ / G ^x- (5.9) 
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Proof, (i) Properties (5.5) and (5.6) follow immediately from the definitions of the three 



tensors and the sets ^ and and are independent of the choice of the reconstruction 
parameters at the interface. 

(ii) Since is assumed to satisfy local force consistency ( |2.7 ), we have 



= (S^UVf) - = |2^l(Sac(l/F;r) - Sa(yF;T)) : drU Wu G %, 



and hence Sac(l/F; •) — ^a{yF] •) is divergence free. According to Lemma 4.4 there exists a 
function ip^'^ G Ni(^)^, which is unique up to a constant shift, such that (5.7) holds. Property 



(5.8) uniquely determines the shift. 
As a matter of fact, it is highly non-trivial whether (5.8 ) can be satisfied, and it is in principle 



possible that the corrections "propagate" into the continuum region [5]. We postpone the 



detailed computations required to prove this to Appendix 6A_ and 6.2[ where we then also give 
a proof of the Lipschitz property (5.9). □ 



5.2. The modified a/c stress. The function ^'^'^(F; •) obtained in Lemma 5.3 provides the 



divergence-free corrector for Sac ~ for homogeneous deformations. We now construct the 
corrector for nonlinear deformations: First, for each / G ^x, f = Tj^ P\T^., we set 

^Av) '■= Wr+y + dT^y). 
We can now define the corrector function for y G as 



J2 r''{ffiyy,mf)Cf+ Yl ^''(2/;^/)C/, 



and the corresponding modified stress function 

Sac(l/;r) := Eac(z/;r)-9V^^^(|/;T)J, 



for T G ^. 



(5.10) 



(5.11) 



We show in Remark 6J^, that ip^'^ is non-trivial, that is, there exists no choice of parameters 
for which ip^'^ = 0, even under purely homogeneous deformations. 

The properties of the modified stress function Sac are summarized in the following lemma. 

Lemma 5.3. Let Sac be defined by ( |5.11[ ), and y G then the following identities hold: 



tUy;T) = i:,{y;T) yre^M 

Sac(l/; T) = S2(y; T) VT G ^Tc; and 
Sac(2/F; •) = Sa(yF; •) VF G M^^l 
Moreover, there exists a constant L^, which depends only on C , such that 



(5.12) 

(5.13) 
(5.14) 
(5.15) 



|Sac(y;T) - Sac(yF;T)| < LacM2p2y| 



^°°(Tn£) 



VT G 



dTy. 



(5.16) 
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Proof. Identity (5.12) follows from (5.2) and the fact that Sac — Sac is divergence-free. 

Identity (5.13) follows from ( 5.5[ ) and the fact that ip^^^y; mj) = for all / G t^^, which im- 
plies that Sac(y; T) = Sac(y; T) = Sa(?/; T) for all T G Similarly, (5.14) follows from (5.6), 
and the fact that ip^^ = ip"^^ in all elements T G 



Fix F G 



1)2x2 



To prove (5.15) we first note that, since ip {y\ 



0, we have tp^^{yf] 



ijj^^(f;»). Using (5.7), we obtain 

Sac(l/F; •) = S 



, F- 



9V'^=(F;.)J = Sa(yF 



We are only left to prove the Lipschitz property (5.16). With F := dry, we have 

|Sac(i/;T)-Sac(i/F;r)| < | Sac T) - Sac (z/f; T) | + 1 9^^^ ; T) - 9^^^^ (|/f; T) | . (5.17) 

From its definition ( 5.3[ ), and the fact that second partial derivatives of V are globally bounded, 
it is clear that Sac satisfies a Lipschitz property of the form 



|Sac(y;7') - Sac(?/F;7')| < Li Afs || /^^Z/ 1| £«>(rn£) 



(5.18) 



where Li depends only on C; see also [9l Lemma 19] for a similar result. (If the reconstruction 
parameters satisfy the one-sidedness condition (3.1), as well as the patch test conditions (2.6), 
(2.7), one may show that Li = SC/Qq.) 

To bound the second term on the right-hand side in (5.17) we invoke the inverse inequality 

2 

^0 



\dr''iy;T)-dr''{yf;T)\ < \r''iy;mf)-r'iy,;mf)\ 



fCT 



where we used the fact that \dCf \ = 2/Qo for all / G If / G then ip" 
f G then tjj^^(»;mf) = il)"^^ m j) and hence, using ( |4.23 ), 



0. If 



If / G J^xi then i/j^'^{y;mf 
employ (5.9) to estimate 



^P''iy;mf)-ij'''{yF;mf)\ = \^^\y;mf)\ < |Af2p'y||£-(rn£)- 

ip'^^{V f]nif) and ip^^{yf;mf) = ip'^^{V]raf). We can therefore 



|V'""(y;m/) - V^"'=(i/f;"^/)| = |^/'"'(Fy;m;) -^/'"=(F;mj)| 



< LacMo F 



/ 



^\<^\\D'y\\nTnc). 



The last inequality can be verified through straightforward geometric arguments. Without 
explicit constants its validity is obvious. 

Combining the two foregoing estimates, we obtain 



\dr\y\T) - dr\yr.T)\ < cmax(Lac, l)M2||D2y||,.o(^nT). 



Combining (5.17), (5.18) and (5.19), yields (5.16). 



(5.19) 

□ 



5.3. Proof of Theorem 5.1 , With the preparations of the foregoing sections it is now easy 
to complete the proof of the main consistency result. Theorem bA_ Again, we drop the depen- 
dence on y whenever possible. We begin by splitting the consistency error into a continuum 



21 



contribution and an interface contribution, 



= J2 1^1 Pac(T) - S,(r)] : + ^ |T| [Sac(r) - Sa(T)] : 9tM 



and estimate E^ and Ej separately. Note also that we used (5.13) to drop the sum over 
elements in the atomistic region. 

Using the fact that Sac = S^ in (5.14), and the stress estimate (4.17), we obtain 



Ec< 5^ \T\\j:l{T)-J:,{T)\\dTu\ 

:(^J2[M,\D''y{x)\ + Ms\DMx)\Ty^\ J2 \T\\dTu\^'^'^'' 



TeSc 
<ci 



%p{C) 



(5.20) 



To estimate Ej, we employ the Lipschitz property (5.16) for Sac and the fact that Sg 

] is also st 

VT G ^, 



under homogeneous deformations (see (5.15)). Using (2.8) it is also straightforward to prove 



dry. 



(5.21) 



Using (5.15), (5.16) and (5.21), we obtain, for any T G Sr, 



|Sac(l/;T) - Sa(i/;r)| < |Sac(i/;T) - Sac(i/F;r)| + |Sa(y;T) - Sa(i/F;T)| 

< L^M2\\D'^y\\eo.(^Tnc) + f^^2||^^?/||f-(Tn£), 
and summing over T G .5^ yields 



Ex< 5^ |T|(La 



^)M2\\D'^y\\(,^^Tnc) \dTu\ 
A i/p' 



< cCxM2\\D'^y\\eP(xc.t)(^Y^ \T\ {druf 



(5.22) 



where Cx depends only on Lac, which depends only on C. 



Combining (5.22) and (5.20) we finally obtain the desired consistency error estimate (5.1). 



This concludes the proof of Theorem 5.1 



6. Appendix: Proof of Lemma 5.2 (ii) 



In this appendix, we provide the remaining details for the proof of Lemma 5.2 (ii). Through- 
out this proof we fix a homogeneous deformation yp, F G M^^^, and drop the argument y = yf 
whenever possible. For example, we will write Sac(T) = Sac(yF;T). 

We begin by computing an expression for Sac — Sa in terms of the parameters Cxj- Equation 
(3.8), in the proof of Lemma 3.2, can be rewritten in the form 



d,v:'{Fsi) = (1 - a,,-i)VF,,-i + c,,jVf, + (1 - a,,+i)^F,,+i. 
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Figure 7. Visualisation of the flat interface analysis in ^6.1 



Recalling also (4.12) and using a,- = a,_i + a,+i, we obtain 



1 ^ 

— 5^ [(1 - a,.^,,_i)rF,,-i + (a,,,,, - l)Vf,j + (1 - a,,,,i+i)VF,,+i] ® a, 
1 ^ 

— J2 ^Fj ® [(1 - j)aj_i + (a^^^, J - l)aj + (1 - j)aj+i] 



1 ^ 



(6.1) 



The explicit evaluation of (6.1 ) for interface elements is carried out separately for fiat interfaces 



and interfaces with corners. 

For triangles not intersecting the interface, Sacd/F;?") — Sa(?/F;T) = 0, hence we need to 
compute the stress errors only for interface elements. 



6.1. Flat interface. Consider the flat interface configuration in Figure |7| According to ( [3lo| ) 
and (3.11) the free parameters are Cj := Cxj (for x G X and j G {2,3,5,6}), and di,i G Z, 
where di = Cxj,^ ^,1 = C^j.^ 4,4, and so forth. We calculate the a/c stress for the elements Ti, T2, 
and collect the results in Table [U 

From Table [1] we can read off the stress differences Sac ^ in the elements Ti,T2: 



Sac(Ti) - Sa(Ti) = {(§ - rf2)VF,l + (c^ - l)Vf,2 + (f " C,)Vf^s} 
+ {{d, - |)Vf,i + (I - C2)V,,2 + (C3 - i)V^F,3} 



0-2 

no 
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Table 1. Table of coefficients of Vfj in (6.1), in interfacial triangles, on fiat interfaces. 



3 

C2 

C3 
2 



a,_i + Cx 



-C3a2 - 



0202 - 
- 0303 

f ia4 



(I - (i2)a2 
(-§ + C2)(a2 - 03) 
(i - C3)(a2 - 03) 
(3 - di)a3 





3 

C3 



3 

C2 
C3 

di 

C5 
2 
3 



C2 
2 



^202 
' 0203 

- |a4 

- 3«5 







C2 
C3 



C2 

C3 
2 

1 
3 

C6 



^1)03 





dia4 — |a5 = (| — rfi)a2 
i"5 ~ 3°6 = (c5 - |)ai 



-dia3 
-C5a4 

-fas + lae - ceOi = (4 - C6)ai 



and 



Sac(T2) - S,(T2) 



= {(rfi-i)VF,i + 

= {(rfi-f)VF,l + 
-{(dl-|)VF,l + 



■2 
3 
■2 
3 
■2 
3 

IC2- 



-C5)l^F,2 + (c6-i)V^F,3} 
-C2)rF,2 + (c3-|)VF,3} 
-C2)I^F,2 + (c3-|)Vf,3} 
- C5)Vf,2 + (C6 - C3)V>,3} 



0.2 

03 
01 



+{(< 

Note that we have provided two alternative representations of Sac(T2) — Sa(T2), since the first 
representation is in general insufficient to construct the corrector. 

Since the atomistic region is a mirror image of the continuum region with respect to the 
interface, we can obtain stress function Sac(l/F; ■) for T3 and T4 from symmetry considerations: 

Sae(T4) - Sa(T4) = {(1 - do)\^F,l + {c, ' l)Vf,2 + (1 - C6)Vf,3} 



+ - 1)V^F,1 + (1 - C5)V>,2 + (C6 - 1)V^F,3} 



0.2 

no 
03 

Uo' 



and 



Sac(r3) - Sa(r3 



= {idi 
= {idi 

-{K 



Uo 
02 

no 



1)Vf,i + (1 - C2)Vf,2 + (C3 - 1)Vf,3} » 
1)Ff,i + (1 - C5)V>,2 + (C6 - 1)Ff,3} ® 
1)Ff,1 + (1 - C5)Ff,2 + (C6 - 1)Vf,3} ® ^ 
-{(C2 - C5)Vf,2 + (C6 - C3)Vf,3} 



Ol 

no- 



From the proof Lemma 4.5| recall that d(f{T)J = —q^clj if / is an edge of T and aj the 
counter-clockwise direction of the edge (relative to T) . We can therefore choose ip^'^ explicitly, 
for example, for / = Ti fl T2: 

(6.2) 



r'{f;mf) := l{{d, - D^F,! + (i - C2)V^F,2 + (c3 - i)V^F,3}. 

For the remaining edges, similar choices can be made, the crucial observation being that the 
terms in neighbouring elements associated with an edge cancel each other out. 

We observe, moreover, that for the triangles Ti and T4, the ai components of the stresses 
vanish, which means that ip^'^{F]mf) = for all / G U <^c- This proves (5.8) in the fiat 
interface case. 
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It remains to prove the Lipschitz bound (5.9). From (6.2) (and the corresponding formulas 
for the remaining edges), it is straightforward to show that ip'^'^ is Lipschitz continous for any 
fixed set of parameters with a Lipschitz constant of th e form LM2, where L can be bounded 
in terms of C. This concludes the proof of Lemma 5.2 (ii) in the flat interface case. 



Remark 6.1 (Correctors are neccessary). From the calculation in this section, it is 
clear that one cannot choose parameters such that T,ac{yF',T) = T,s.{y^]T) for all T ^ ,'7' and 
for all potentials V E Y. For example, if Eac(?/F;T2) = Sa(?/F; ^2) for all V, then di = 2/3, 
whereas if Sac(2/F;23) = Sa(|/F;T3), then di = 1. This demonstrates that the divergence-free 
corrector fields are in fact necessary, and that it is impossible in our current framework to 
construct an a/c method where Sacd/F;^) = ^a.{yf',T) holds for all T G ,'7,F E M?^^, and 

V eY. □ 



6.2. General interface. We now turn to the proof of (5.7)- (5.9) for interface configurations 
with corners. Consider the corner configuration displayed in Figure |8| which is concave from the 
point of view of the atomistic region. The reconstruction coefficients found in Proposition 3^ 
are displayed in the figure as well. Recall that the reconstructions of bonds into the atomistic 
or continuum regions are now uniquely determined, while the bonds lying at the interfaces 
(parameters a and b) are still free. 

Using (6.1), and defining a'j := ttj/flo, the stress errors Sac — Sa in the elements Ti 



can again be computed explicitly: 



Sac(Ti) 


- Sa(Ti 


Sac(T2) 


- Sa(T2 


Sac(T3) 


- Sa(T3 


Sac (^4) 


- Sa(T4 


Sac (^5) 


- Sa(T5 


Sac(^6) 


- Sa(T6 



) = (|Vf,3 - |Vf,2) 8) a[ + (a - |)Vf,i ® 4 - (a - f )Vf,i ® a'. 



3; 



(fe- i)V>,3®a'i, 

- (6 - |)\/f,3 ® a[ + {b- f )Ff,3 ® 4 + (|Vf,i - 3 VF,2 

(|Vf,2 - |V>,i) 8) 4 + [(1 - a)\/F,i + {b- l)rF,3] ® 4 

- (6 — 1)Vf,3 ® «i5 S'^d 

(|Vf,2 - |V>,3) ® a'l + [{a - l)V^F,i + (1 - &)Vf,3] ® 4 

- (a - 1)Ff,i (8)03. 



3' 
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Figure 9. All possible corner configurations (up to translation, rotation and reflection) 



Following the argument in ^6.1 , we can check again that the associated edge contributions from 



neighbouring elements cancel, and hence we can explicitly construct the corrector function 



Note that Sac(T2) has no a'^ component and Sac(T'5) has no a'^ component, which implies (5.8). 

For a corner that is convex from the point of view of the atomistic region, the result follows by 
symmetry (interchanging the coefficients 1 and |). The Lipschitz bound (5.9) can be obtained 



from the above formulas, under the assumption that the reconstruction coefficients a, b are 
bounded above by C. 

Finally, we have to convince ourselves that our above argument applies to all possible in- 
terface geometries. In Figure |9] we present an exhaustive list, up to translations, rotations 
and reflections, of local interface geometries. (Recall our geometric requirements formulated 



in Assumption 3.5 ) By inspecting the calculation of the stress differences Sac — Sa for the 
case presented in Figure [8} one observes that the formulas are local, and do not depend on 
the extended geometry of the interface. We note, however, that this only holds due to the 



separation Assumption |3.5[ The subsequent construction of the corrector now follow of course 
verbatim. 



This concludes the proof of Lemma 5.2 (ii) in the general interface case. 



7. Conclusion 

We have shown for a 2D model problem that it is possible to construct patch test consistent 
a/c coupling method for multi-body potentials, in interface geometries with corners, using a 
new variant of the geometry reconstruction technique introduced in [3l [17] , which we labelled 
the GR-AC method. Moreover, we have proven a quasi-optimal consistency error estimate for 
the GR-AC method(s) we constructed. 

We see this work as a first step towards a general theory of GR-AC method(s). Our goal is to 
show eventually that the free parameters in the method can always (that is, in any dimension, 
for any interface geometry) be determined so as to satisfy the energy and force consistency 
conditions, and that the resulting GR-AC method(s) will have the same consistency properties 
that we establish in the present case. 

An important issue that we have left entirely open in the present work is the stability of the 
GR-AC method: Under which conditions on the reconstruction parameters does the GR-AC 
method have sharp stability properties as discussed in [2]? This issue is the topic of ongoing 
research. 
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